package gov.fnal.mcginnis.threetone;

import java.io.File;

/**
 * @author mcginnis
 *
 */
public abstract class SimplexOptimize 
{
	int	MAX_NO_OPT_VARIABLES = 100;

	int		m_inum_total_variables;
	int		m_inum_opt_variables;
	double	m_fsim_limit;
	int		m_iopt_limit_type;
	COptimizeVariable[]	m_variable = null;
	int[]	m_ilookup = null;
	File	m_FwriteAmoebaText;
	
	class COptimizeVariable
	{
		String	cname;
		double  fvalue;
		double	fmin;
		double	fmax;
		double	fstep;
		double	finit_value;
		double	imask;
	}
	abstract double	fChiSquare();
	abstract void	setupOptimization();
	int iAmoeba(double[][] p, double[] y, int ndim, double ftol,
					int[] nfunk)
	{
		int		i;
		int		ihi;
		int		ilo;
		int		inhi;
		int		j;
		int		mpts=ndim+1;
		double	rtol;
		double	sum;
		double	ysave;
		double	ytry;
		double[]	psum = null;
		double	ftime = 0.0;
		double  dswap;
		double	TINY = 1.0e-10;
	   
		psum= new double[ndim + 1];
		nfunk[0] =0;
		for (j=1;j<=ndim;j++)
		{
			for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];
			psum[j]=sum;
		}
		for (;;) 
		{
			ilo=1;
// converted weird NR code
			if( y[1] > y[2] )
			{
				ihi = 1;
				inhi = 2;
			}
			else
			{
				ihi = 2;
				inhi = 1;
			}
			for (i=1;i<=mpts;i++) 
			{
				if (y[i] <= y[ilo]) ilo=i;
				if (y[i] > y[ihi]) 
				{
					inhi=ihi;
					ihi=i;
				} 
				else if (y[i] > y[inhi] && i != ihi) inhi=i;
			}
			if ( (nfunk[0] > ((int) m_fsim_limit)) 
				&& ( m_iopt_limit_type == 1 ) )
			{
	/* NMAX exceeded */
				dswap = y[1];
				y[1] = y[ilo];
				y[ilo] = dswap;
				for (i=1;i<=ndim;i++) 
				{
					dswap = p[1][i];
					p[1][i] = p[ilo][i];
					p[ilo][i] = dswap;
				}
				psum = null;
				return -1;
			}
			if ( (ftime > m_fsim_limit) && ( m_iopt_limit_type == 2 ) )
			{
	/* Time exceeded */
				dswap = y[1];
				y[1] = y[ilo];
				y[ilo] = dswap;
				for (i=1;i<=ndim;i++) 
				{
					dswap = p[1][i];
					p[1][i] = p[ilo][i];
					p[ilo][i] = dswap;
				}
				psum = null;
				return -2;
			}
			if ( (y[ilo] < m_fsim_limit) && ( m_iopt_limit_type == 3 ) )
			{
	/* Chi min exceeded */
				dswap = y[1];
				y[1] = y[ilo];
				y[ilo] = dswap;
				for (i=1;i<=ndim;i++) 
				{
					dswap = p[1][i];
					p[1][i] = p[ilo][i];
					p[ilo][i] = dswap;
				}
				psum = null;
				return -3;
			}
			rtol=  (2.0*Math.abs(y[ihi]-y[ilo])/(Math.abs(y[ihi])+Math.abs(y[ilo])+TINY));
			if (rtol < ftol) 
			{
				dswap = y[1];
				y[1] = y[ilo];
				y[ilo] = dswap;
				for (i=1;i<=ndim;i++) 
				{
					dswap = p[1][i];
					p[1][i] = p[ilo][i];
					p[ilo][i] = dswap;
				}
				break;
			}
			nfunk[0] += 2;
			ytry=amotry(p,y,psum,ndim,ihi,-1.0);
			if (ytry < 0.0)
			{
				psum = null;
				return -100;
			}
			if (ytry <= y[ilo])
			{
				ytry=amotry(p,y,psum,ndim,ihi,2.0);
				if (ytry < 0.0)
				{
					psum = null;
					return -100;
				}
			}
			else if (ytry >= y[inhi]) 
			{
				ysave=y[ihi];
				ytry=amotry(p,y,psum,ndim,ihi,0.5);
				if (ytry < 0.0)
				{
					psum = null;
					return -100;
				}
				if (ytry >= ysave) 
				{
					for (i=1;i<=mpts;i++) 
					{
						if (i != ilo) 
						{
							for (j=1;j<=ndim;j++)
								p[i][j]=psum[j]= 0.5 * (p[i][j]+p[ilo][j]);
		
							vSetupTrial(psum);
							y[i] = fChiSquare();
							if (y[i] < 0.0)
							{
								psum = null;
								return -100;
							}
						}
					}
					nfunk[0] += ndim;
					for (j=1;j<=ndim;j++)
					{
						for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];
						psum[j]=sum;
					}
				}
			} else --(nfunk[0]);
		}
		psum = null;
		return 0;
	}
	double amotry(double[][] p, double[] y, double[] psum, int ndim,
				   int ihi, double fac)
	{
		int		j;
		double	fac1;
		double	fac2;
		double	ytry;
		double[] ptry = null;
	   
		ptry= new double[ndim + 1];
		fac1=( 1.0 - fac)/ndim;
		fac2=fac1-fac;
		for (j=1;j<=ndim;j++) ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
		
		vSetupTrial(ptry);
		ytry = fChiSquare();
		if (ytry < 0.0)
		{
			ptry = null;
			return ytry;
		}
		if (ytry < y[ihi]) 
		{
			y[ihi]=ytry;
			for (j=1;j<=ndim;j++) 
			{
				psum[j] += ptry[j]-p[ihi][j];
				p[ihi][j]=ptry[j];
			}
		}
		ptry = null;
		return ytry;
	}
	int Optimize(double ftol)
	{

		int		isolution;
		int		iparam;
		int[]	infunc = new int[2];
		int		istatus;
		double[]	fvary_param = null;
		double[]	fchi_square = null;
		double[][]	fsolutions = null;
		double	fchiba;

		fvary_param = new double[m_inum_opt_variables + 1];
		fchi_square = new double[m_inum_opt_variables + 1  + 1];
		fsolutions = new double[m_inum_opt_variables + 1  + 1][m_inum_opt_variables  + 1];

		fchiba = 0.0;
		for ( iparam = 1; iparam <= m_inum_opt_variables; ++iparam )
		{
			fvary_param[iparam] = 0.0;
		}

		vSetupTrial(fvary_param);
		fchiba  = fChiSquare();
		if (fchiba < 0.0)
		{
			fvary_param = null;
			fchi_square = null;
			fsolutions = null;
			return -1;
		}
		for (	isolution = 1; 
				isolution <= m_inum_opt_variables + 1; 
				++isolution )
		{
			for ( iparam = 1; iparam <= m_inum_opt_variables; ++iparam )
			{
				fvary_param[iparam] = 0.0;
				fsolutions[isolution][iparam] = 0.0;
			}
			if ( isolution > 1 ) 
			{
				fvary_param[isolution - 1] = 1.0;
				fsolutions[isolution][isolution - 1] = 1.0;
			}
			vSetupTrial(fvary_param);
			fchi_square[isolution] = fChiSquare();
			if (fchi_square[isolution] < 0.0)
			{
				fvary_param = null;
				fchi_square = null;
				fsolutions = null;
				return -1;
			}
		}
		istatus = iAmoeba(fsolutions,fchi_square,
						m_inum_opt_variables,ftol,infunc);
/*		
		if ( istatus == -1 )
		{
			sprintf(cmessage,"Exceeded maximum number of iterations %d",infunc);
			vWriteAmoebaText(cmessage);
		}
		if ( istatus == -2 )
		{
			sprintf(cmessage,"Exceeded Time limit of %f",m_fsim_limit);
			vWriteAmoebaText(cmessage);
			sprintf(cmessage,"Number of iterations %d",infunc);
			vWriteAmoebaText(cmessage);
		}
		if ( istatus == -3 )
		{
			sprintf(cmessage,"Reached Chi Square of %f",m_fsim_limit);
			vWriteAmoebaText(cmessage);
			sprintf(cmessage,"Number of iterations %d",infunc);
			vWriteAmoebaText(cmessage);
		}
		if ( istatus == 0 )
		{
			sprintf(cmessage,"Function minimized after %d iterations",infunc);
			vWriteAmoebaText(cmessage);
		}
		if ( istatus == -100 )
		{
			sprintf(cmessage,"Minimization interrupted");
			vWriteAmoebaText(cmessage);
		}
*/
		vSetupTrial(fsolutions[1]);
		fchiba = fChiSquare();
		
		fvary_param = null;
		fchi_square = null;
		fsolutions = null;
		
/*		if (fchiba > 0.0)
		{
			sprintf(cmessage,"Ending Chi Square = %f",fchiba);
			vWriteAmoebaText(cmessage);
			return 0;
		}
		else
		{
			sprintf(cmessage,"Minimization interrupted");
			vWriteAmoebaText(cmessage);
			return -1;
		}
*/
		return istatus;
	}
	void vSetupTrial(double x[])
	{
		int		ii;

		for (ii = 0; ii < m_inum_opt_variables; ++ii )
		{
			m_variable[m_ilookup[ii]].fvalue 
				= x[ii + 1] 
				* m_variable[m_ilookup[ii]].fstep
				+ m_variable[m_ilookup[ii]].finit_value;
			if (m_variable[m_ilookup[ii]].fvalue 
				> m_variable[m_ilookup[ii]].fmax)
			{
				m_variable[m_ilookup[ii]].fvalue 
					= m_variable[m_ilookup[ii]].fmax;
			}
			if (m_variable[m_ilookup[ii]].fvalue 
				< m_variable[m_ilookup[ii]].fmin)
			{
				m_variable[m_ilookup[ii]].fvalue 
					= m_variable[m_ilookup[ii]].fmin;
			}
		}
		return;
	}
	void vWriteAmoebaText(String cmessage)
	{
/*		if (m_FwriteAmoebaText == NULL) return;
		fprintf(m_FwriteAmoebaText,"%s\n",cmessage);
		fflush(m_FwriteAmoebaText);
*/	}

}
